setwd("~/Dropbox/research/losingtouch/replication") # Remember to set this to your working directory

library(tidyverse)
library(haven)
library(skimr)
library(estimatr)
library(huxtable)
library(ggeffects)
library(scales)
library(viridis)
library(modelsummary)
library(fixest)

# Define colors
show_col(viridis(9))
vircol1 <- viridis(9)[1]
vircol2 <- viridis(9)[2]
vircol3 <- viridis(9)[3]
vircol4 <- viridis(9)[4]
vircol5 <- viridis(9)[5]
vircol6 <- viridis(9)[6]
vircol7 <- viridis(9)[7]

dir()

options(scipen=10000)

sd <- read_dta("experiment-data.dta")

sdrlong1 <- sd %>% 
  dplyr::select(id,starts_with("topic"),polpref) %>% 
  pivot_longer(topic_a:topic_b) %>% 
  mutate(option=str_sub(name,-1)) %>% 
  transmute(id,option,topic=value)

sdrlong2 <- sd %>% 
  dplyr::select(id,starts_with("easy"),polpref,gender,profile_age2,profile_education) %>% 
  pivot_longer(easy_a:easy_b) %>% 
  mutate(option=str_sub(name,-1)) %>% 
  transmute(id,option,polpref,easier=value,
            gender=factor(gender),agegrp=factor(profile_age2),edu=factor(profile_education))

sdrlong <- sdrlong2 %>% 
  left_join(sdrlong1,by=c("id","option")) %>% 
  mutate(polpref2=case_when(option=="b" & polpref!=97 ~ polpref,
                            option=="a" & polpref!=97 ~ (polpref-10)*-1,
                            TRUE ~ NA_real_))


easylm1 <- lm_robust(polpref2~easier+topic,data=sdrlong,clusters = id)
summary(easylm1)

#print regression table for full model
easylm0 <- feols(polpref2~easier,data=sdrlong,cluster = "id")
easylm1 <- feols(polpref2~easier | topic,data=sdrlong,cluster = "id")
easylm2 <- feols(polpref2~easier | topic,data=sdrlong,cluster = "id",subset = sdrlong$topic!="Schools")

modlist <- list(`Model 1` = easylm0,`Model 2` = easylm1, `Model 3` = easylm2)

modelsummary(modlist,
             output = "latex",
             fmt = 3,
             booktabs = TRUE,
             estimate = "{estimate}{stars}",
             statistic = "({std.error})",
             coef_map = c("easier" = "Easier snippet",
                          "(Intercept)" = "Intercept"),
             gof_omit = c("AIC|BIC|R2|RMSE"))

easylm2 <- lm_robust(polpref2~easier+topic,data=sdrlong,clusters = id,subset=topic!="Schools")
summary(easylm2)

# average preference for easier and harder snippets
elm2pr90 <- ggpredict(easylm2,terms = "easier",ci.lvl = .90) %>% 
  as_tibble() %>% 
  transmute(conf.low90=conf.low,conf.high90=conf.high)

elm2pr <- ggpredict(easylm2,terms = "easier") %>% 
  as_tibble() %>% 
  bind_cols(elm2pr90) %>% 
  mutate(cond=fct_inorder(c("Complex","Simple")))

elm2dvsd <- sdrlong %>% filter(topic!="Schools") %>% pull(polpref2) %>% sd(.,na.rm=T)
elm2eff <- tidy(easylm2)[2,2]
elm2eff/elm2dvsd #[1] 0.1144834

ggplot(elm2pr) +
  aes(x=cond,y=predicted) +
  geom_errorbar(aes(ymin=conf.low,ymax=conf.high),width=0,color=vircol3,size=.5) +
  geom_errorbar(aes(ymin=conf.low90,ymax=conf.high90),width=0,color=vircol3,size=1) +
  geom_point(color=vircol3) +
  theme_bw() +
  labs(x="Proposal style",y="Predicted preference (0-10)") +
  theme(axis.text.x = element_text(color="black",size=10))

ggsave("figures/expplot.pdf",width=6,height = 4)
ggsave("figures/expplot_pres.pdf",width=7,height = 3)
ggsave("figures/expplot_tall.pdf",width=5,height = 6)

